function F = cal_chg_bartik(param, const, eqm)
%% Choose optimal quality after Bartik shock
        
    % Define PiB
    DB_check = ((eqm.D_H).^(param.beta_v/(param.beta_v-1)) + (param.e^param.sigma*eqm.D_F*param.shock).^(param.beta_v/(param.beta_v-1))).^((param.beta_v-1)/param.beta_v);
    D1_check = ((eqm.D_H).^(param.beta_v/(param.beta_v-1)) + (param.e^param.sigma*eqm.D_F).^(param.beta_v/(param.beta_v-1))).^((param.beta_v-1)/param.beta_v);
    PiB = (DB_check./D1_check).^const.gamma.*eqm.Pi1; 
        
    % Export decision
    ExportCutoffQB = (const.k0*const.lnzQ-log(param.sigma*const.gamma)+log(kron(ones(param.Nf,1),PiB)-kron(ones(param.Nf,1),eqm.Pi0))-param.mu_E)/param.sigma_E;
    ExportProbQB = normcdf(ExportCutoffQB);         
    
    % Grid search
    [qB_index, ~] = grid_search(param, const, eqm.Pi0, PiB, ExportCutoffQB, ExportProbQB); 

      
%% Calculate average wage response

    % Quality level
    q_level = const.Qgrid(eqm.q_index)'; 
    qB_level = const.Qgrid(qB_index)';
    
    % Firm wage
    wageO = const.wage_grid(eqm.q_index)';
%     wageB = const.wage_grid(qB_index)';
    
        % used in earlier estimation
        pp_lnwage = interp1(log(const.Qgrid(eqm.nq~=0)), log(const.wage_grid(eqm.nq~=0)), 'spline', 'pp');
        hat_lnwage = @(lnq) ppval(pp_lnwage, lnq);
        hat_wage = @(q) exp(hat_lnwage(log(q)));   
        wageB = hat_wage(qB_level);  
    
    
    % Average quality and wage response  
    weight = eqm.density_exp;
    avg_quality_response = sum(qB_level./q_level.*weight)/sum(weight);    
    avg_wage_response = sum(wageB./wageO.*weight)/sum(weight);
 
 
%% Changes in average wage of suppliers   
  
    avg_wgt_lnwageS_q = eqm.theta_m.*eqm.c.^(param.sigma-1)./eqm.V.*sum(const.phi_y.*const.phi_v.*kron(ones(param.Q,1), (log(param.w)+const.lnwage_grid).*eqm.Psig),2)'; 
    avg_wgt_lnwageS = avg_wgt_lnwageS_q(eqm.q_index)';
    avg_wgt_lnwageSB = avg_wgt_lnwageS_q(qB_index)';
    avg_wageS_response = sum((avg_wgt_lnwageSB - avg_wgt_lnwageS).*weight)./sum(weight);

    
%% Export intensity

    rB = (eqm.D_H/param.f_vH).^(1/(param.beta_v-1))./((eqm.D_H/param.f_vH).^(1/(param.beta_v-1))+(param.e^param.sigma.*eqm.D_F.*param.shock/param.f_vF).^(1/(param.beta_v-1)));
    DB = rB.*eqm.D_H + (1-rB).*param.e^param.sigma.*eqm.D_F.*param.shock;     
    export_intensity = 1-rB.*eqm.D_H./DB;        


%% Return

    bartik.PiB                  = PiB;
    bartik.qB_index             = qB_index;
    bartik.avg_quality_response = avg_quality_response;
    bartik.avg_wage_response    = avg_wage_response;
    bartik.avg_wageS_response   = avg_wageS_response;
    bartik.export_intensity     = export_intensity;
        
    F = bartik;

end